Distribution of particles and bubbles in turbulence at small Stokes number 
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The inertia of particles driven by the turbulent flow of the surrounding fluid makes them prefer 
certain regions of the flow. The heavy particles lag behind the flow and tend to accumulate in the 
regions with less vorticity, while the light particles do the opposite. As a result of the long-time 
evolution, the particles distribute over a multi-fractal attractor in space. We consider this distribu- 
tion using our recent results on the steady states of chaotic dynamics. We describe the preferential 
concentration analytically and derive the correlation functions of density and the fractal dimensions 
of the attractor. The results are obtained for real turbulence and are testable experimentally. 
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Recently the problem of distribution of inertial par- 
ticles in turbulence received much attention of the re- 
searchers [H-H2I]. This is largely thanks to the break- 
through in the theoretical understanding of the La- 
grangian motion of particles in the flow that occurred 
lately [13j. While the understanding of the behavior of 
particles that have negligible inertia and follow the flow 
is quite complete by now El, the understanding of the 
behavior of inertial particles is still insufficient. This is 
while the subject has an extremely wide range of applica- 
tions: the flows of fluids are typically turbulent and often 
laden with external particles. Theoretical advancement 
was made mainly for the case of small Stokes number, 
where the inertia is weak and the particles " almost" fol- 
low the flow. Even in this limit of small Stokes numbers, 
the particles' distribution is highly non-trivial. Particles' 
deviations from the surrounding flow accumulate with 
time, bringing particles to a strange attractor in space. 
This attractor is multi-fractal and the only theoretical 
result obtained so far for the real turbulent flow was 
the derivation of the correlation codimension 0|. Here 
a result obtained for real turbulence is a result obtained 
without modeling turbulence and expressed in terms of 
the (unknown) statistical properties of turbulence. Since 
the statistics of turbulence is largely unknown [l4| , then 
to obtain such a result one needs to make universal pre- 
dictions on particles' behavior in the flow independent of 
the details of the statistics of that flow. 

In this Letter we provide the complete description of 
the distribution of particles in real turbulence at small 
Stokes numbers, describing both the correlation of the 
particles' density with the surrounding flow and the 
statistics of the singular density on the attractor. We give 
a number of predictions that are testable experimentally. 

The idea that particles' inertia leads to inhomogeneous 
spatial distribution dates back to the seminal paper by 
Maxey 0. It was observed that due to inertia heavy 
particles are pushed out of the vortices and hence they 
will not distribute uniformly in the flow, like the inertia- 
less particles. However the quantitative description of 
the correlations between the locations of particles and 
of vortices stayed unaddressed. Note the distribution of 
vorticity in turbulence is random and dynamical, while 



the distribution of particles reflects its cumulative effect 
over time. There is a residual correlation that we describe 
by an integral relation holding in the steady state. 

We find the spectrum of fractal dimensions of the at- 
tractor. We show that while the correlation dimension is 
different from the dimension of space, the fractal or sim- 
ilarity dimension [l5| is equal to the space dimension. In 
contrast the information dimension is different from the 
spatial dimension and it equals the Kaplan- Yorke dimen- 
sion. In turn, the correlation codimension equals twice 
the Kaplan- Yorke codimension which constitutes a pre- 
diction allowing direct testing in the laboratory. 

The analysis is based on the recent finding of a univer- 
sal description for the steady state density of the weakly 
compressible dynamical systems [16| . The particles' mo- 
tion, though governed by Newton's law, admits an effec- 
tive description in terms of a velocity field in space. Iner- 
tia is described by a small compressible correction to the 
incompressible velocity of the background turbulent flow. 
This correction leads to a small disbalance of trajecto- 
ries going in and out of space regions, which accumulates 
over a long time to a big effect. Thus compressibility is 
a singular perturbation which treatment was performed 
in [l6|. For a mixing incompressible velocity the evo- 
lution of a small volume of particles makes it dense in 
space. The volume's coarse-graining over an arbitrarily 
small scale covers all the available space, which volume 
is assumed finite. When a small compressible compo- 
nent is added to the velocity, the coarse-graining of the 
evolved volume over an arbitrarily small scale does not 
cover the whole space any longer. However the coarse- 
graining over a small but finite scale, that tends to zero 
with compressibility, already covers the whole volume. 

The analysis assumes the single-particle approximation 
where one neglects the interaction between the particles 
and their back reaction on the flow. We consider a small 
spherical particle with the radius a and the material den- 
sity p p suspended in a fluid with the density p and the 
kinematic viscosity v. The fluid flow u(t,x) is assumed 
to be incompressible. The Newton law governing the 
evolution of the particle's position x(t) and the particle's 
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velocity v(t) is assumed to have the form 

dv d , , S1 v — u\t, x(t)] 

where 7 = 3p/(p + 2p p ) and r = a 2 /(3^7) is the Stokes 
time. Thus we assume that all the forces besides the 
added mass and the drag can be neglected @, ■ With 
no loss we set the total volume and the mass equal to 
one, so the spatial average of the particles' density n 
obeys (n) = 1. We set j3 = 7 — 1 so the particle's velocity 
relative to the flow w(t) = v(t) — u[t, x(t)] obeys 



dw w „ d , 



(2) 



The parameter (3 — 2(p — p p )/(p + 2p p ) changes from —1 
for heavy particles to 2 for light ones. After transients 



w(t) = (i 



exp 



t-t' 



dt 



7 u[t',x(t')]dt'. (3) 



We assume r is much smaller than the smallest time- 
scale of turbulence, which is the viscous time-scale t v , 
so the Stokes number St = T/t v <C 1. Then we can 
substitute the derivative in the integrand by its value at 
t! = t so v{t) « u + p [d t u + (u ■ V)u] with p = j3r — 
2a 2 (p—p p )/(9vp). Thus at St -C 1 the particle's velocity 
is determined uniquely by its position x(t) in space and 
one can introduce the particle's velocity field v(t, x) 

x(t) = v [t, x(t)] , v = u + p [dtu + (u ■ V)ii] . (4) 

In the zero inertia limit St —> the particles follow the 
incompressible mixing flow of turbulence x = u[t,x(t)] 
and in the steady state they are uniformly distributed in 
space, so their steady state density n s equals one. This 
behavior is characteristic of small dye particles. However, 
at a small but finite St, the small correction v — u gives 
the particles' velocity field a finite compressibility [2] 



W = V • V = -fl(j) ^ 0, <p = lo 2 - s 2 



(5) 



so the constant is no longer a solution to the continu- 
ity equation 9 t n + V ■ (nv) = 0. Above s 2 = SijSij and 
to 2 = aijdij, where Sij is the symmetric (strain) and ay 
is the antisymmetric (vorticity) parts of velocity gradient 
djUi = Sij + a.ij. The field <f>(x) is positive in the regions 
dominated by vorticity and negative in the regions dom- 
inated by the strain, and it will be called below the indi- 
cator, indicating whether x is in a vortex. It follows from 
the Navier-Stokes equations that <f> equals the Laplacian 
of the turbulent pressure, 4> = V 2 p. Eq. ([5]) shows that 
heavy particles (3 < are repelled from vortices (here 
and below "vortex" is used qualitatively), while the light 
ones /3 > are attracted. This is the generalization of 
the familiar fact that a heavy particle in a centrifuge is 
pushed out to the boundary. Turbulence can be con- 
sidered as a dynamically changing spatial distribution of 
vorticity, so heavy particles tend to accumulate on the 



boundaries between the vortices, cf. [H, 0,1111: forming a 
singular density supported on these boundaries. This ac- 
cumulation however is insignificant during the life-time of 
a single vortex and the ultimate singular distribution of 
particles in space n s forms from the long-time combined 
action of many uncorrelated vortices. Still one expects a 
residual correlation between the distributions of vortic- 
ity and particles, to find which we consider the steady 
state density n s . One expects n s can be obtained by let- 
ting an arbitrary initial condition n in the remote past 
n(t = — T) = no evolve for infinite time, T — > 00, accord- 
ing to the continuity equation. Starting from the uniform 
initial distribution we obtain the steady state density 



n s (x)—\im n(T), n(T) = exp 

T— s-oo 



)[t, q(t, x)]dt 



(•0) 



provided a condition of decay of correlations [16j . that 
should hold in our case, is met. For the cross-correlation 
of density and vorticity F(x) = (cj)(0)n s (x)) we find 



F(x)=( 0(0,0) exp 



p, <f>[t,q(t,x)]dt 



(7) 



d a In ( exp 



a<t>(0,0)+(M / cf>[t,q(t,x)}dt 



\a=Q, 



where we used the conservation of the mean de nsity 
(n(t)} = const. Applying the cumulant expansion [171 ] . 
taking derivative of the series and setting a = 0wc find 

F(x) = p / dt((4(0, 0)#, q(t, x)])) e + 0(St 2 ).(8) 
J — 00 

The above formula is the same as one would obtain by 
expanding the exponent in Eq. |(7J| and keeping the lowest 
order term in r, with one important difference. In Eq. ([S]) 
one has the second order cumulant or dispersion, that one 
would not get by the series expansion of the exponent, cf. 
Q. This difference is essential as without the cumulant 
the integral in Eq. © diverges: (w[t,q(t,x)]) = (1 — 
/3)r(0[t, q(t, x)]} is equal to a non-zero sum of Lyapunov 
exponents, see below. To leading order in St one can 
substitute q(t,x) in Eq. ([5]) by X(t,x) 

d t X(t,x) =u[t,X(t,x)}, X(0,x) = x. (9) 

where X(t, x) are Lagrangian trajectories of u. One finds 

{mn s {x))=pf dtM0,0)#,X(t,aj)]), (10) 

J — 00 

where we can already omit the cumulant since by 
incompressibility (w[t,X(t,x)]) = (w(t,x)), while 
/ w(t,x)dx = J V • v(t,x) = by the boundary con- 
ditions. Since there is no degeneracy, the non-negative 
spectrum of the Laplacian of pressure in the Lagrangian 
frame (j)\t, X(t, x)] is strictly positive at zero frequency 



25(0) 



M0,0)#,X(t,0)]>>0. (11) 
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The single-point correlation (<pn s ) = J <fr(x)n s (x)dx 
equals [iE(0)/2 and it gives the integral of <j> where each 
region weighted by the number of particles in it 



J[uj 2 {x)-s 2 (x)] n s {x)dx = a 2 (p- p p ){$vp)- 1 E(0),( 



12) 



where we used the definitions of <p and p. For heavy 
particles, p p > p, the answer is negative giving a measure 
of the extent to which the particles favor regions with 
negative cf>. For light particles, p p < p the answer is 
positive measuring their favoring of vortices. The above 
integral steady state relation holds at any t. 

The quantity E(0) appeared first in [J], where it was 
shown to determine the correlation dimension of the par- 
ticles' attractor in space. This quantity is increased by 
the intermittency of turbulence and it can be estimated 
as t~ 3 f(Re) where /(-Re) is a growing function of the 
Reynolds number Re that grows as a power d, Q . We 
show E(0) determines all the fractal dimensions. 

We observe that Eq. Q is a weakly dissipative dy- 
namical system, defined as the dynamics for which the 
potential part of v is much smaller than the solenoidal 
one. The statistics of the steady state density of such 
systems was shown recently to allow for a complete and 
universal description [16j . The application to our case 
gives the following results. The motion of particles in 
space is chaotic and is characterized by the Lyapunov 
exponents 19]. To the lowest order in St the exponents 
are equal to the Lyapunov exponents Xi of the turbulent 
flow u. However, the value of the sum of the Lyapunov 
exponents A^~ that determines the logarithmic rate of 
growth of the volumes forward in time, 



;A+(x) = lim t' 1 IndetV jqi{t,x) 



(13) 



is zero for it, so the leading order approximation demands 
the account of the correction v — u. This is also the 
case of the sum of the Lyapunov exponents ^2 X[ of the 
backward-in-time flow that determines the density [l6j 



lim T~ x ]nn(0,x,T) = X^(x). 



(14) 



the Kaplan- Yorke codimension Cky of the attractor. At 
weak compressibility the definition [20( of Cky reduces 
[IH to Cky — I ^3 which gives to the leading order 
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KY 



P 2 E(0) = 2a\p-p p ) 2 E(0) 
2|A 3 | 8I1/VIA3I ' 



(15) 



where the third Lyapunov exponent A3 determines the 
rate of exponential separation of X(t,x) back in time. 
We have |A 3 | ~ t~ l and C KY ~ P 2 St 2 f(Re), cf. 0,0]. 

The probability for two particles to be at the dis- 
tance x is described by the pair-correlation function 
(n s (0)n s (x)). Substituting for n s the expression from 
Eq. ([5]) and using the cumulant expansion [16j one finds 
(n s (0)n s (x)) = exp [/x 2 g(a;)] where the structure func- 
tion g(x) depends only on the statistics of turbulence 



g(x)= f <ft 1( fta<0[ti,X(ti,O)]^[t2,X 

J —OO 



(*2,x)])(16) 



The above is valid if the higher order terms in the cu- 
mulant expansion are negligible [16]. The Kolmogorov 
theory (KT) estimate would give the validity condition 
St <C 1, while the account of intermittency changes the 
condition to h{Re)St <C 1 where h(Re) is expected to be 
a slowly growing function of Re. cf. [||. The function 
g(x) has a universal behavior [16[ at small x that gives 



(n s (0)n s (x)) = (r)/x) 



2Ck 1 



x <C 77, 



(17) 



where r\ ~ {vt^ 1 ! 2 is the Kolmogorov scale of turbulence 
[H . Thus the correlation codimension equals ICky ■ 

Remarkably, the structure function determines all the 
correlation functions of n s . Generalization of the calcu- 
lation of {n s (Q)n s (x)) gives the log-normal statistics [lH ] 



(n s (x 1 )n s (x2)-n s (x N )) =cxp 



. (18) 



The density n s does not have physical meaning and we 
consider the coarse-grained density n;, 



For turbulence J2^t is expected to be the same for all mi{x) = j n s (x')dx', ni(x) = 3m;(a;)/(47rZ 3 ), 

x with the possible exception of a set of points with zero J\x'-x\<i 



x with the possible exception of a set of points with zero 
volume. The results of [l|| HH give 



(w(0,0)w[t,X(t,0)] 



p 2 E(0) 



where mi(x) is the mass in a small ball. The limits I — > 
and St do not commute ((nf) - (77 / 'l) 2C ' KY ') 



Thus for all initial points, with the possible exception of 
a set of points with zero volume, the infinitesimal vol- 
umes decay to zero in the limit of infinite evolution time, 
while the steady state density n s is zero except for a set 
of points with zero volume. Due to conservation of mass 
J ndx, we conclude that n s has S— function type singu- 
larities on its support. This support is the strange at- 
tractor - the multifractal set in space that is approached 
by the particles' trajectories at large times. We now find 



lim lim (n> ) 



1. 



lim lim(n ; ) 

St->0i->0 



(19) 



For any St > the fluctuations of ni are large for a suf- 
ficiently small I. On the other hand, by lim.st^o(nf) = 1 
one sees that for any fixed / > the fluctuations of ni are 
small for a sufficiently small St. The coarse-grained den- 
sity is uniform over scales which minimal value vanishes 
with St. Thus turbulence effect on the particles depends 
on the observer's resolution I: at 2Cky^(v/^) ^ 1 se S" 
regation holds, while 2Cky In (?7/0 < 1 ■ mixing. This 



4 



is how mixing works effectively for particles on a multi- 
fractal. Segregation may also bring physical effects [4| . 

At St <C 1 there is a scale L <C 77 over which the 
density is almost uniform. We note that m;(i = 0, cc) 
is equal to the mass contained in the preimage of the 
ball time t ago, which is an ellipsoid around q(—t,x) 
with the largest axis growing as Zexp[|A3f|]. At t* = 
A3I -1 ln(L/Z) the ellipsoid has the scale over which the 
density is uniform, so the mass contained in it is just its 
volume AttI exp [-/_, w[t',q(t',x)]dt']/3 and we find 



ni(x) — exp 





-/< / <j,[t',q{t',x)]dt' 

-\\ 3 \-i\n{L/t) 



(20) 



see [16| for details. The smallness of fi brings the ex- 
pected conclusion that the statistics of ni is log-normal 



(v/i) 



Cky p(p-l) 



(21) 



cf. Eq. (fT5)) . The spectrum of the fractal dimen- 
sions D(a) = lim;^.o ln(m" _1 n s )/[(a — l)lnZ] in- 
volves the average with n Sl rather than the spatial 
average 0, [T3|. To find it consider (n" _1 n s ) = 
limT^oo(exp[— a f_ t uj[t,q(t,r)]dt — J_^"ui[t,q(t,r)]dt]). 
Due to St < 1 the contribution of time-intervals 
with length t v is negligible and we may substitute 
the upper limit in the last integral by — — t n 
which allows to perform independent averaging 
(exp[— a J_ t u[t, q(t, r)]dt — tv uj[t,q(t,r)]dt]) « 

(exp[— a f_ t oj[t,q(t,r)]dt]}{exp[—f_^* tv uj[t,q(t,r)]dt]}. 
However the last average is equal to one by the 
conservation of mean density, so = (nf) and 



D(a) = 3 - Ckycl. 



(22) 



Our results generalize to the two-dimensional case, where 
they can be compared with [6J]. Working out the small 



compressibility limit reproduces our answer. Returning 
to the three-dimensional case, we observe that the fractal 
dimensions are close to 3 (we do not consider a> 1) The 
fractal dimension of the attractor D(0) coincides with the 
space dimension 3, which is somewhat counter- intuitive 
since the volume of the attractor is zero. The information 
dimension -D(l) is equal to the Kaplan- Yorke dimension. 

At St <C 1 density inhomogeneities are absent in the 
inertial range [l4|. Then Eq. (fl7)) is a complete descrip- 
tion. In contrast, at St ~ 1 the inertial range inhomo- 
geneities are important Q and Eq. (fl"6"]) . extended to 
hold asymptotically at St ~ 1, gives a unique access 
to the inhomogeneities. In KT g{x) depends only on 
x and the mean energy dissipation e so (n s (0)n s (x)) — 
exp [C/i 2 e 2 / 3 a;~ 4 / 3 ] . This prediction describes correctly 
the model of v decorrelated in time [3|, [3, |8| . However, for 
turbulence, simulations @ show \n{n s (0)n s (x)) oc x~ 10 / 3 
at moderate Re where KT is expected to work. Noting 
g{x) ~ T^d^(\p(x)— p(0)] 2 ), where t x is the relevant time- 
scale, we suggest the difference has the same origin as the 
deviations of the pressure scaling from KT [2l| . 

A central result is the analytic description of the prefer- 
ential concentration by Eq. (fT2"|) . A single number E(0) 
completely characterizes the influence of turbulence on 
the log-normal statistics of density at r <C rj. Log- 
normality arises because the steady state density is the 
cumulative result of the creation of inhomogeneities by 
many uncorrelated vortices, each of which creates but 
weak inhomogeneity. The fractal structure at scale I 
forms relatively fast - within the characteristic time-scale 
A3 1 — 1 \n(i]/l). The predictions are testable. 
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